Surfing on the Edge: 
Chaos vs. Near-IntegrabiUty in the System of Jovian Planets 



Wayne B. Hayes 

Computer Science Department University of California, Irvine 
Irvine, California 92697-3435 

wayneSics . uci . edu 
ABSTRACT 

We demonstrate that the system of Jovian planets 
(Sun+Jupiter+Saturn+Uranus+Neptune), integrated for 200 milhon years 
as an isolated 5-body system using many sets of initial conditions all within 
the uncertainty bounds of their currently known positions, can display both 
chaos and near-integrability. The conclusion is consistent across four different 
integrators, including several comparisons against integrations utilizing quadru- 
ple precision. We demonstrate that the Wisdom- Holman sym plectic map using 



simple symplectic correctors as implemented in Mercury 6.2 (IChambersI Il999l ) 
gives a reliable characterization of the existence of chaos for a particular initial 
condition only with timesteps less than about 10 days, corresponding to about 
400 steps per orbit. We also integrate the canonical DE405 initial condition out 
to 5 Gy, and show that it has a Lyapunov Time of 200-400 My, opening the 
remote possibility of accurate prediction of the Jovian planetary positions for 5 
Gy. 

Subject headings: outer solar system, dynamics, chaos 



1. Introduction 

Both the man of science and the man of art live always at the edge 
of mystery, surrounded by it. Both, as a measure of their creation, have 
always had to do with the harmonization of what is new with what is familiar, 
with the balance between novelty and synthesis, with the struggle to make 
partial order in total chaos.... This cannot be an easy life. 

— J. Robert Oppenheimer 
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When one speaks of the stabihty of our Solar System, one must carefully define the 
meaning of "stable" . We say that the Solar System is practically stable if, barring interlopers, 
the known planets suffer no close encounters between themselves or the Sun, over the main- 
sequence lifetime of the Sun. In a practically stable Solar System, the orbital eccentricities, 
inclinations, and semi-major axes of all the planets remain within some bounded region, not 
too far from their present values. In this sense, work by many au thors over the past 15 years 



has all but proven that the Solar S ystem is practically s table (jLaskar 



Laskai]ll997l : llto and Tanikawalbond ). Good reviews exist (iLissauerl 



1999 



994 



Laskai 



Lecar et al. 



1996 



200ll) 



and we will not discuss it further in this paper. A second, more formal definition involves 
the question of whether the Solar System is chaotic or not. In a chaotic system, nearby 
solutions tend to diverge from each other exponentially with time, although in a weakly 
chaotic system such as the Solar System, the exponential divergence can be preceded by 
an initial period of polynomial divergence. Let d{t) be the distance between two solutions, 
with d{0) being their initial separation. Then d{t) increases approximately as (i(O)e^* in a 
chaotic system, where A is the Lyapunov exponent. The inverse of the Lyapunov exponent, 
1/A, is called the Lyapunov time, and measures how long it takes two nearby solutions to 
diverge by a factor of e. A system that is not chaotic is called integrable or regular, and has 
a Lyapunov exponent of zero. A practical consequence of being chaotic is that small changes 
become exponentially magnified, so that uncertainties in the current positions of the planets 
are magnified exponentially with time. Even though the Solar System is practically stable, a 
positive Lyapunov exponent means that uncertainties in the current positions of the planets 
are magnified to the point that we cannot predict the precise positions of the planets in their 
orbits after a few (or at most a few tens of) Lyapunov times. 

KAM theory tells us that essentially all Hamiltonian systems which are not integrable 
are chaotic. An initial condition (IC) not lying precisely on a KAM torus will eventu- 
ally admit chaos, but with a t i me s c ale that depends critically on the IC. Symp lectic 
integrators ( Channell and Scovell Il990l : I Wisdom and HolmanI Il99ll : ISanz-Sernal Il992l ) have 
many nice properties when used for long-term integrations of Hamiltonian systems, such 
as conservation of phase-space volume, and bounded energy error. However, the validity 
of symplectically-integrated numerical solutions also depends critically upon the integration 
tin ie step h, with the longevity of the solu tion's validity scaling as e"/'* for some constant 
a (IBenettin and Giorgillil Il994j : iReichI Il999l ). For linear problems, the dependence is even 
stronger and manifests itself as a bif urcation i n the L yapunov exponent , going discontinu- 
ously from zero to a non- zero value (jLessnickl (119961 ). iNewman and Led (120051 ) — but see 
Ranch and HolmanI (119991 )). Since the solar system is not integrable, and experiences unpre- 
dictable small perturbations, it cannot lie permanently on a KAM torus, and is thus chaotic. 
The operative question is the time scale of the chaos. To compute the time scale accurately. 
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we must guarantee that the measured time scale is not an artifact of the integration method. 



What is the Lyapunov time of the Solar System? ISussman and Wisdom! (119881 ) first 
demonstrated that the motion of Pluto is chaotic with a Lyapunov time of about 2 million 
years, corroborated over a longer integration later by iKinoshita and Nakail (119961 ). iLaskar 
( 1l989l ) performed an averaged integration of the 8 major planets (excluding Pluto) and 
found that the Lyapunov ti me was about 5 million years, with the divergence dominated 
by that of the inner planets. iLaskan (119901 ) beli eved that secular resonance s are the cause 
of the chaos in the inner Solar System (but see iMurray and HolmanI (120011 ) ) , although he 
did not believ e the system of Jovian plane ts was affected by the chaos displayed by the 
inner planets. ISussman and Wisdom] (119921 ) performed a full (non-averaged) integration of 
the entire Solar System and confirmed Laskar's 5 million year Lyapunov time, and further 
found that the system of Jovian planets by itself had a Lyapunov time of between 7 and 
20 million years, although their measurement of the Lyapunov time displayed a disturbing 
dependence on the timestep of the integration. This dependence was later discovered to 
be due to symplectic integratio n schemes effectively integrating a slightly di fferent set of 
ICs; the effect can be corrected (ISaha and Tremaindll992l : IWisdom et al.lll996l ). although it 
decreases with decreasing timestep. 

Since there are no two-body resonances am ongst the Jovian pla n ets, t he cause of the 
chaos between them was not understood until iMurrav and HolmanI (Il999l) iden tified the 
cause as being the overlap of three-body resonances. iMurray and HolmanI (119991 ) also per- 
formed Lyapunov time measurements on a large set of Outer Solar Systems differing only in 
the initial semi-major axis of Uranus. They found that their three-body resonance theory 
correctly predicted which regions of ICs were chaotic, and which were not, at least over the 
200-million-year integration timespan they used. For the "ac t ual" S olar System, they found 
that the Lyapunov time was about 10 million years. iGuzzol (12005! ) went on to corroborate 
the three-body resonance theory by performing a large suite of integrations, numerically 
detecting a large web of three-body resonances in the outer Solar System. 



Murray and HolmanI (jl999l ) noted that the widths Aa/a of the individual resonant zones 



was of order 3 x 10^^, so that changes in the ICs of that order can lead to regular motion. They 
note, however, that "the uncertainties in the ICs, and those introduced by our numerical 
model, are comfortably smaller than the width of the individual r esonances, so [the outer] 
Solar System is almost certainly chaotic." Given that iGuzzol ( l2005l ) has also detected many 
three-body resonances consistent with Murry -|- Holman's theory, it would seem at first 
glance that chaos in the outer Solar System is a fact. 



However, the conclusion that the isolated outer Solar System is chaotic cannot be taken 
for granted. For example, it is known that symplectic integration with too-large a timestep 
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can inject chaos into an integrable system (IHerbst and Ablowitzlll989l : iNewman et al.ll2000l ). 
Although most authors verify their primary results by performing "checking" integrations 
with smaller timesteps, the checking integrations are not always performed for the full du- 
ration of the main integrations. Thi s, combined with the fact that longer sy mplectic in- 
tegrations require shorter timesteps (IBenettin and Giorgillil 1 19941 : iReichI Il999l ) means that 
one cannot assume that a timestep good enough, for example, for a 100-million-year inte- 
gration is also good enough for a 200-million-year integration. There is currently no known 
method for analytically choosing a short-enough timestep a priori, and so the only method of 
verifying the reliability of an integration is to re- perform the ent i re inte gration with shorter- 
and-shorter timesteps until the results conver ge. iNewman et al.l (I20001) us ed this method to 
demonstrate that, for a given set of ICs, the IWisdom and Holmanl (jl99ll ) symplectic map- 
ping with a 400 day timestep (about 11 steps per orbit, a commonly-quoted timestep) admits 
chaos, but that the results converge to regularity for any timestep less than about 100 days. 
However, many authors who find chaos have also performed reasonable convergence tests, 
demonstrating that the chaos does not always disappear at convergence. 

There exists compelling evidence for the absense of chaos in the outer Solar System. 
Laskar and others noted that when the entire solar system is integrated, the inner Solar 
System manifests chaos on a 5-million-year timescale, but the outer Solar System appears 
regular in these integrations. Although Laskar's approximate theory can overlook some 
causes of chaos, there also exist full-scale integrations that indicate the absense of chaos. 
Grazier et al. (1999) and Newman et al. (2000) utilized a Stormer integrator whose per-step 
numerical errors were bounded by the double-precision machine epsilon (about 2 x 10"^^), 
as long as the orbital eccentricity is less than 0.5 and 1000 steps or more are taken per 
orbit. Furthermore, their integration method took care to ensure that the roundoff error was 
unbiased. Note that, except for the possibility of having the same property with a larger 
timestep, this is practically as good an integration as is possible using double-precision. 
Furthermore, such an integration is symplectic by default, since it is essentially exact in 
double precision. Using this method, they perforr ned an integr a tion o f the Jovian planets 
lasting over 800 million years, and found no chaos. IVaradi et al.l (120031 ) performed a 207My 
integration of the entire Solar System, including even the effects of the Moon, and placed a 
lower bound of 30My on the Lyapunov time of the system of Jovian planets. 

We are thus left with the disturbing fact that, utilizing "best practices" of numerical 
integration, some investigators integrate the system of Jovian planets and find chaos, while 
others do not. 



In this paper, we demonstrate that this apparent dilemma has a simple solution. Namely, 
that the boundary, in phase space, between chaos and near-integrability is finer than pre- 



- 5 - 



viously recognized. In particular, the c urrent observat i onal uncertainty in the po sitions of 
the outer planets is a few parts in 10^ (IStandishlll998l : iMorrison and Evanslll998l ). Within 
that observational uncertainty, we find that some ICs lead to chaos while others do not. So, 
for example, drawing 7-digit ICs from the same ephemeris at different times, one finds some 
solutions that are chaotic, and some that are not. Thus, different researchers who draw 
their initial coniditions from the same ephemeris at different times can find vastly different 
Lyapunov timescales. 



2. Methods 



Wit h the exception of the two sets of initial condition s (ICs) we have received from other 
auth ors (IMurrav and HolmanI 1 19991 : I Grazier et al.l |2005| ) and the set included in Mercury 
6.2 (I Chambers Il999l ). all ICs used in this paper are drawn at various epochs from DE405 
(jStandishlll998l ). which is the latest planetary ephemeris publicly available from JPL. It has 
stated uncertainty for the positions and masses of the outer planets of a few parts in 10^. 
To ensure that our integration agrees over the short term with DE405, we verified in several 
cases that we can integrate between different sets of DE405 ICs, separated by as much as 
100 years, while maintaining at least 7 digits of agreement with DE405. 

We integrate the system of Jovian planets using only Newtonian gravity. The inner 
planets are accounted for by adding their masses to the Sun and perturbing the Sun's position 
and linear momentum to equal that of the Sun+Mercucy+Venus+Earth+Mars system. This 
ensures that the resonances between the out er planets is shifted by an a mount that is second 
order in this mass ratio, roughly 3 x 10~^^ (IMurray and Holmanlll999l ). which is far smaller 
than the uncertainty in the outer-planet positions. We assume constant masses for all objects 
and ignore ma ny effects which are probably relevant over a 200My timescale (see for example 



Laskai ( 



( iLaskar 



199 



We account for solar mass loss at a rate of m/m f» 10 ^ per million years 



19991 : lNoerdlingerll2005l ). but note that we observe no noticable difference if we keep 



the solar mass constant. 

To reduce the possibihty that our results are dependent on the integration scheme, 
we used three different numerical integration methods to verify many of our results in this 
paper. First, we used the M ercury 6.2 package (jChamberslll999l ). with the Wisdom-Holman 
( IWisdom and HolmanI Il99ll ) symplectic mapping option (called MVS in the input files). We 
used stepsizes varying from 2 days to 400 days. Second, we used the NBI package, which 
contains a 14th-order Cowell-Stormer method with modifications by the UCLA research 
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group led by William Newman ( Grazier et al. 1995 . 1996 : Varadi et al. 20031 )1^ NBI has been 



shown to have relative truncation errors below the double precision machine epsilon (about 
2 X 10~^^) when more than 1000 steps per orbit are used and the orbital eccentricity is less 
than 0.5. More precisely, if the largest component of the phase-space vector of the solution at 
time t has absolute value M, then the local errors per step of each of the components are all 
less than 2 x lO'^^M. Note that this means that the component-by-component relative error 
can be significantly greater than the machine precision for components of the solution that are 
significantly less than M, but all components have errors relative to M which are smaller than 
the machine precision. Furthermore, the authors of NBI have gone through great pains to 
ensure that the roundoff error is unbiased. We used a 4-day timestep for all NBI integrations, 
which gives more than 1000 timesteps per Jupiter orbit. We have verified the above "exact to 
double precision" property by comparison against quadruple precision integrations described 
below. Note that for the above-defined timestep, NBI gives practically the best integration 
possible using only double precision, and that such an integration is symplectic by default 
since it essentially provides a solution which is exact in double precision. In our 200 million 
year integrations, NBI always had relative energy errors and angular momentum errors of 
less than 2 x 10^^^, with an average of about 2 x 10^^^. 



Our third integrator was the Taylor I.4 package (iJorba and Zoul 120051 ). Taylor I.4 is 
a recent and impressive advance in integration technology. It is a general-purpose, off-the- 
shelf integrator which utilizes automatic differentiation to compute arbitrary order Taylor 
series expansions of the right-hand-side of the ODE. Taylor I.4 automatically adjusts the 
order and stepsize at each integration step in an effort to minimize truncation error, and 
utilizes Horner's rule in the evaluation of the Taylor series to minimize roundoff error. As 
the authors note, integration accuracy is gained more efficiently by increasing the order of 
the integration than by decreasing the timestep, since the accuracy increases exponentially 
with the order but only polynomially in the timestep. Although Taylor I.4 allows the user 
to specify a constant order and timestep, we chose to allow it to use variable order and 
timestep while providing it with a requested relative error tolerance equal to 1/1000 of the 
machine precision, in order to produce solutions which were exact to within roundoff error. 
We found that Taylor I.4 typically used about 27th order with about a 220 day timestep. 
The fact that the solution is exact to machine precision over such a long timestep guarantees 
that accumulated roundoff is by far the smallest in the Taylor I.4 integrations. Furthermore, 
Taylor integrators are extremely stable when applied to non-stiff problems, with the radius 
of convergence increasing linearly with integration order, and in our case the timestep is 



^ NBI is available at Ihttp : //astrobiology . ucla . edu/~varadi/NBI/NBI . html or by searching the web 
for "NBI Varadi". 
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well within the radius of convergence (IBarrio et aLll2005l ). Finally, Taylor I.4 allows the 
user to specify the machine arithmetic to use, including software arithmetics. Out-of-the- 
box, Taylor I.4 supports the use of IEEE 754 double precision (64 bit representation with 
a 53 bit mantissa), Intel extended precision (80 bit representation with a 64-bit mantissa, 
giving a machine precision of about 10~^^, ac cessible as l ong double when using GCC on an 
Intel machine), the DoubleDouble datatype (jBriggslll996l ) which provides software quadruple 
precision in C++, and the GNU Multiple Precision Library, which allows arbitrary precision 
floating point numbers in C++. Most of our integations using Taylor I.4 used Intel extended 
precision, which is almost as fast as double precision and gives about 19 decimal digits 
of accuracy. Over our 200-million-year integrations using Intel extended precision, Taylor 
typically had relative energy errors of less than 8 x 10~^^; the worst relative energy error 
observed in any of our integrations was 2 x 10"^'^. Integrations began with the Solar System's 
barycentre at the origin with zero velocity. After 200 million years the barycentre drifted a 
maximum of 3 x 10~^° AU, while the z component of the angular momentum was always 
conserved to a relative accuracy better than 3 x 10~^^. We also performed a suite of quadruple 
precision integrations, in which energy and angular momentum were each conserved to at 
least 26 signiflcant digits over 200 million years. 
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Brouwer's Law: Taylor 1 .4's error is roundoff-dominated at most stringent tolerance 
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Fig. 1. — Taylor I.4 satisfies Brouwer's Law: when local tolerance is set above machine 
precision of 10~^^, phase-space error grows as due to biased truncation errors (upper 
three curves). But with tolerance set well below the machine precision, its result is exact 
{i.e., correctly rounded) to machine precision, so that phase-space error grows only as t^'^ 
(lowest curve). 
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Now we analyze the error growth as a function of time for our non-symplectic integrators, 
apphed to a non-chaotic IC. One consequence of having a solution whose numerical error is 
dominated by unbiased roundoff (ie., exact to machine precision) is that when integrating a 
non-chaotic system, the total phase-space error grows polynomially as t^'^. If the local error is 
biased (either by biased roundoff, or due to truncation errors in the integration scheme), then 
the total phase-space error grows as t^. (This is assuming that the integrator is not inherently 
symplectic, as is the case with bo th NBI and T a ylor l .j.) This is known as Brouwer's Law 



( lBrouwerlll937l ). As noted above, iGrazier et al.l (l2005l ) have demonstrated that the error in 



NBI's integration is dominated by unbiased roundoff when using 1000 timesteps per orbit. We 
tested Taylor i.^ in Intel extended precision for similar properties by integrating the system 
of Jovian planets for 200 million years using various integration tolerances up to and beyond 
the machine precision of 10~^^, and compared these integrations to a Taylor I.4 integration 
that used quadruple precision. In Figured], we plot the phase-space separation between the 
quadruple precision integration and several integrations using Intel extended precision. We 
see that when the relative integration tolerance is set above the machine precision, the error 
grows as t"^ and is therefore truncation dominated. But when the tolerance is set to 10^^^ 
(about a factor of 1000 below the machine precision), the error grows as t^'^, and is therefore 
dominated by unbiased roundoff. This is consistent with Taylor I.4 producing results that 
are exact in Intel extended precision when given a local relative error tolerance of 10~^^, just 
as NBI produces exact results in double precision when used with 1000 or more timesteps 
per orbit. We see that after 200 million years, the errors in the positions of the planets are 
of order 10~^ AU (in the case of non-chaotic ICs). This translates into a phase error, or 
equivalently an observational erro r, of substanti a lly les s than one arc-minute. A comparison 



of our Figure [T] with Figure 2 of I Grazier et al.l (119991 ) at their 200 million year mark also 



demonstrates that Taylor I.4 using Intel extended precision provides about 3 extra digits 
of precision over NBI using double precision, as would be expected when comparing 19-digit 
and 16-digit integrations. As we show later (Figure [7]), the differences in orbital elements 
are even smaller. 

For comparison, we briefly mention the run-times of the integration algorithms. All 
timings are for a 2.8Ghz Pentium 4 processor. As we shall see later, the largest timestep for 
which the Mercury ^j.^ integrations agreed with the others was 8 days; a 16-day timestep was 
almost as good. Thus we compare the Wisdom- Holman mapping with timesteps of 16 and 8 
days to NBI with a 4-day timestep, and Taylor-l.^- In addition to the Taylor-l.4 extended 
precision (long double) timings used in this paper, we present its timings for standard IEEE 
754 double precision, for comparison to the other double precision integrations. Table [T] 
presents the results. We found that the total runtime was linearly proportional to the inverse 
of the timestep, as would be expected. We note several observations. First, Taylor-l.4 is 
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not competitive in terms of efficiency. However, note that Taylor-1.4 is a "proof-of-concept" 
software package for general-purpose integration of any system of ODEs, and it currently 
generates code that can be quite inefficient. A carefully hand-coded Taylor series integrator 
for Solar System integrations is far more efficient, and can be competitive with the above 
codes (Carles Simo 2007, personal communication). Second, Wisdom-Holman with an 8- 
day timestep is the fastest case among the integrations we tested that showed complete 
convergence. Third, Wisdom-Holman with a 4-day timestep (51 hours, not shown in the 
table) is slower than NBI with a 4-day timestep. Finally, although we did not test NBI 
for convergence at timesteps less stringent than 4 days, it is possible that NBI maintains 
the lead in being more efficient than Wisdom-Holman, if it also shows convergence at larger 
timesteps. 
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integrator 


WH dt=16d 


WH dt=8d 


NBI 


Taylor- 1.4 double 


Taylor-1.4 long double 


time (hours) 


13 


26 


40 


100 


150 



Table 1: Run times for integrating the 5-body system for 200 million years on a 2.8Ghz 
Pentium 4, for Mercury 6.2 with timesteps of 16 and 8 days, for NBI, and for Taylor-1.4- in 
double and long double precision. 
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We will not directly measure the Lyapunov time in this paper, since it is difficult to 
create an objective measure of the Lyapunov time over a finite time interval. Formally, 
the Lyapunov time is only defined over an infinite time interval. In practical terms, the 
divergence is almost always polynomial for some nontrivial duration before the exponential 
divergence emerges, and it is difficult to pinpoint the change-over objectively. However, 
when plotting the distance between two nearby trajectories as a function of time, it is usually 
evident by visual inspection whether or not exponential divergence has occured by the end of 
the simulation. Thus we will plot the actual divergence between two numerical trajectories 
initially differing by perturbing the position of Uranus by 10~^^ AU (about 1.5mm) in the z 
direction. We will call these pairs of trajectories "siblings." In the cases where we see only 
polynomial divergence between siblings over a 200 million year integration, we will abuse 
terminology and call these systems "regular", "ncar-intcgrablc" , or "non-chaotic", although 
formally all wc have shown is that the Lyapunov time is longer than can be detected in a 
200-million-year integration. 
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Survey of chaos while varying Uranus semi-major axis, Quad Precision 
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Fig. 2. — Murray + Holman (1999) in quadruple precision. 



Survey of chaos while varying Uranus semi-major axis, Intel 80-bit Precision 
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Fig. 3.— Murray + Holman (1999) in long double. 
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As a prelude to our main results, we first crudely reproduce the results of iMurray and Holman 



( 119991 ) ■ Murray + Holman performed a survey in which the semi- major axis of Uranus ajj 
was varied around its actual value, plotting the measured Lyapunov time as a function of 
ajj. Broadly, they found that below ajj ~ 19.18, the Lyapunov time was substantially less 
than 10'' years, sometimes as small as 10^ years; in the range au ~ 19.20 19.22, the Lya- 
punov time fluctuated between about 1-10 million years, interspersed with cases in which 
the motion was regular; near au = 19.24 it was uniformly regular; near au = 19.26 there was 
a tightly-packed region of both chaotic and regular orbits; and near au = 19.28 the motion 
was again uniformly regular. We have reproduced this survey using quadruple precision in- 
tegrations, using a very course grid in au due to computing constraints. (On a 2.8GHz Intel 
Pentium 4, a quadruple precision integration using Taylor I.4 proceeds roughly at lOOMy 
per month of CPU time.) Figure [2] displays the distance between siblings, for various values 
of au- Although we were not able to perform the survey using a finer grid in au due to the 
computational expense, we see that in broad outline we obtain results similar to Murray -|- 
Holman. In Figure [3l we repeat the same survey using Intel extended precision. We see 
that the results are qualitatively identical, demonstrating that 19 digits of precision (which 
requires about 1/20 of the CPU time of quadruple precision) gives qualitatively identical, 
and quantitatively very similar, results. Henceforth in this paper, all Taylor I.4 integrations 
are performed using Intel extended precision, with the relative local error tolerance set to 

10-22. 



3. Results 
3.1. Corroborating previous results 

As an e arly step, the author obtain ed from Murray + Holman their initial conditions 



(ICs) used in [Murray and Holman! (119991 ). and verified using accurate integrati ons that their 



system was chaotic. The author then obtained from P. Sharp the ICs used by iGrazier et al. 



(I2OO5I ). and verified that their system was regular, at least over a 200My timespan. The 
author then spent significant time eliminating several possible reasons for the discrepancy, 
such as incorrect ICs, incorrect deletion of the inner planets, etc. The author also integrated 
from one set of ICs to the other (having along the way to account for the fact that they were 
in different co-ordinate systems), and finally found that the systems agreed with each other 
to a few parts in 10^ when integrated to the same epoch. After some time it became apparent 
that neither group of authors had made any obvious errors. The conclusion seemed to be 
that both systems nominally represented the outer Solar System to within observational 
error, but one was chaotic and one was not. 
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3.2. Ephemeris Initial Conditions Drawn at Different Times 

It is reasonable to make the assumption that either our Solar System is chaotic, or it 
is not. It cannot be both. This appears to be the assumption that most practitioners make 
when measuring "the" Lyapunov time of our Solar System. Although this is probably a 
reasonable assumption, it does not follow that all initial conditions (ICs) drawn from an 



ephem eris are equivalent. The most recent ephemeris published by JPL is DE405 (jStandish 



19981 ). It is based upon hundreds of thousands of observations, all of which have finite error. 



Thus, the ephemeris does not represent the exact Solar System. For the outer pl anets, the 



best match of DE405 to the obse rvations yields residual errors of a few parts in 10^ (jStandish 



19981 : iMorrison and Evanslll998l ). The masses of the planets are also known to only a few 
parts in 10^. It is possible that within these error bounds there exist different solutions with 
different Lyapunov times. In particular, it may be possible that some solutions display chaos 
on a 200-million-year timescale, while others do not. 
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Fig. 4. — Distance in AU between "sibling" trajectories as a function of time for 21 sets of 
initial conditions drawn from DE405 at 30-day intervals. Each plot represents one initial condition. 
Siblings are integrated with the Wisdom-Holman symplectic map with timesteps: 50 (solid line), 
100 (long dash), 200 (short dash), and 400 days (dotted). For any given plot, note disagreement 
across timesteps, and that the switch between chaos and regularity is not always monotonic with 
timestep. 
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To test this hypothesis, we drew 21 sets of ICs from DE405, at 30-day intervals starting 
at Juhan date 2448235 (9.5 Dec, 1990). We represent each of our 21 sets of ICs using a 3-digit 
numeral, 000, 030, 060, . . ., 570, 600, representing the number of days after JD2448235 at 
which the intitial condition is drawn. We chose a 2-year total interval across which to draw 
our ICs in order to to ensure a reasonable sample of inner-planet positions before deleting 
the inner planets (a Martian year is about 2 Earth years). We drew the initial conditions by 
taking the output of the program testeph.f, which is included with the DE405 ephemeris. 
The output of testeph.f is rounded to 7 digits, since no more digits are justifiable!! As 
described above, we augment the mass of the Sun with that of the inner planets and augment 
the Sun's position and linear momentum to match that of the inner Solar System. For each of 
the 21 ICs, we generate a "sibling" IC which is offset by 10"^^ AU (1.5mm). We then integrate 
both of them for 200 million years, and measure the distance between them at 1-million- 
year intervals. To ensure that our results are not integrator- or timestep-dependent, we 
perform each integration in several ways: firs t, with the Wisdom- Holman symplectic mapping 
included with Mercury 6.2 Jchamberslll"999[ ) with timesteps of 400, 200, 100, 50, 32, 16, 8, 4, 
and 2 days; second, with NBI with a timestep of 4 days; and third, with Taylor I.4 in Intel 
extended precision with a relative error tolerance of 10~^^. Results for the Wisdom-Holman 
integrations with larger timesteps are displayed in Figure ID After staring at these figures 
for awhile, several key observations become apparent. First, when comparing across the 21 
sets of initial conditions, there is remarkable disagreement about whether or not the outer 
Solar System displays exponential divergence. This will be discussed further below, using 
more accurate integrations. Looking at an individual graph, but across timesteps, we note 
that there are very few cases (e.g. t = 150, 330, 420) in which there is universal agreement 
between all four timesteps that divergence is polynomial. There are also only a few cases 
that universally agree that the divergence is exponential. In most cases, there is substantial 
disagreement across timesteps whether a particular initial condition admits chaos. This is 
consistent with the observation of Newman et al. (2000). However, in contrast to Newman 
et al. (2000), we note that the "switch" from chaos to non-chaos is not always monotonic in 
timestep. For example, for system 000, timesteps of 400 and 50 days admit chaos, while the 
"in-between" timesteps of 200 and 100 days do not. For system 180, timesteps of 400 and 
200 days display non-chaos, while for 100 and 50 day timesteps, chaos is apparent; this is 
precisely opposite to what one would expect if large timesteps were injecting chaos into the 
system. As we shall see later in the paper (Table [2]), there appears also to be no observable 
correllation between the timestep and the percentage of ICs that admit chaos. Thus we 



^Although this is probably not the best way to uniformly sample initial conditions from within the error 
ball representing the error in the observations, it doe.s represent a reasonable way to reproduce how users of 
DE405, taking initial conditions from DE405 to 7 digits, get their samples. 
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hypothesize that the discrepancy across timesteps is due more to perturbations in the IC s 
and our use of only low-order symplectic correctors (IWisdom et al.lll996l : IChamberd Il999l ) . 
than it is due to unreliable integration at large timesteps. Corroborating this hypothesis 
would require us to re-perform these ex periments using higher-order symplectic correctors 
or "warmup" (jSaha and Tremainelll992l ). which is a possible direction for future research. 
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Fig. 5. — Similar to previous Figure, but more accurate integrations. NBI: solid lines; Wisdom- 
Holman: dashed lines with dt =8 days (long dash), 4 days (short dash); Dotted lines {Taylor-1.4) 
are shifted below others because Taylor-1.4 is more accurate. For any given IC, there is good 
agreement between the shapes of the curves, indicating agreement on the existence (or lack) of 
chaos. (Divergence saturates at 100 AU.) 
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Fig. 6. — Curves identical to the previous figure, but plotted on a log-log scale, which depicts 
polynomial divergence as a straight line. 
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Figures [5] and [6] plot the divergence between sibling trajectories for all 21 ICs, as inte- 
grated by more accurate integrations showing convergence: NBI, the Wisdom-Holman map- 
pings with timesteps of 8 and 4 days, and Taylor 1.4- Wisdom-Holman with 2 day timesteps 
agreed with these curves, but are omitted to reduce clutter; Wisdom-Holman with a timestep 
of 16 days also showed good agreement in all but two cases. Both figures are identical except 
that Figure [5] uses a log-linear scale, while Figure [6] uses a log-log scale; different features 
are visible using the different scales. After staring at these graphs for awhile, at least two 
observations present themselves. First, if one looks at the system corresponding to any sin- 
gle IC, there is usually good agreement between the integrations as to the future divergence 
between the sibling trajectories of that particular case. This demonstrates that convergence 
has occured and makes it unlikely that the results are integrator dependent. Second, looking 
across cases, it appears that the future of the outer Solar System over the next 200 million 
years is quite uncertain, varying from nearly integrable, to chaotic with a Lyapunov time of 
order 10 million years or less. This is quite a startlingly diverse array of possible outcomes, 
considering that the ICs for these systems are all drawn from the same ephemeris, all less 
than 2 years apart, and presumably differing from each other by only a few parts in 10''. 
The author has, in fact, verified that several of the ICs, when all integrated up to the same 
epoch in the vicinity of 1991, agree with each other to a few parts in 10''. 
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Fig. 7. — The final difference at tlie 200My mark between siblings, in semi-major axis and 
eccentricity, for the 21 cases depicted in the previous Figures. The integrator was NBI. 
Eccentricity difference is the Euclidean distance between siblings in the 4-dimensional space 
consisting of the four orbital eccentricities. Semi-major axis difference is the Euclidean 
distance between siblings in the 4-dimensional space consisting of the four Aa/a values. 
Some values have been moved slightly for textual clarity, but in no case by more than the 
width or height of a character. 
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Figures [5] and [6] demonstrate that in the chaotic cases, the sibhngs can have their 
respective planets on opposite sides of the Solar System after 200My. However, Figure [7] 
demonstrates that changes in the orbital elements are much less drastic, demonstrating that 
the Solar System is practically stable over a 200My timescale even when it is chaotic. 



3.3. Explicitly Perturbed Initial Conditions 



Following [Murray and HolmanI (119991 ) , we performed several surveys in which we per- 
turbed the semi-major axis au of Uranus from its current value, but kept all other initial 
conditions (ICs) constant. We used all three previously mentioned integrators; Wisdom- 
Holman with timesteps of 4 and 8 days; NBI; an d Taylor- 1. A. Th e IC was the default one 
from the file big. in included with Mercury 6.2 (jChamberd Il999l ) . which according to the 
documentation is from JD2451000.5. The inner planets were deleted, with their mass and 
momentum augmenting the Sun's as described elsewhere in this paper. We completed sur- 
veys in which au was changed in steps of 2 x 10^^ for k = 6,7, 8, which corresponds to 
Aau/o-u in steps of lO"''^^^'*. We went 10 steps in each direction for each value of k. For 
each step, we generated a "sibling" IC by randomly perturbing the positions of all planets 
by an amount bounded by 10~^^ Aujfl We then integrated both for 200My, and plotted the 
distance between them as a function of time. For k = 7,8, there was no significant difference 
between any of the integrations. That is, all siblings at all steps had virtually identical 
divergences when changing au in 10 steps of 2 x 10"''^^'^^ AU, for a total change of 2 x 10~^. 
This corresponds to Aau/au stepped by 10"''^*'^^, for a total change in Aau/au of 10"'' in 
each direction. However, for k = 6, some of the steps showed chaos while others did not. 
The change was not monotonic: over the 21 steps (10 in either direction plus the "baseline" 
case), there were three "switches" between chaos and stability. 



Note that, for no good reason, this is different from the perturbations used to generate sibhngs in the 
rest of the paper. 
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Two systems with delta a_U/a_U = 1e-7, one chaotic, one not 
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time in millions of years 

Fig. 8. — Two systems, both using ICs within the error bounds of the best known positions 
for the outer planets. Both have identical initial conditions except for the semi-major axis 
au of Uranus, which differs between the two systems by 10~^ in ^au/au- One admits chaos, 
while the other does not. The "upper" six curves (all starting with sibling distances near 
10~^) are all double-precision integrations, two using Wisdom-Holman with timesteps of 4 
and 8 days, and one using NBI. Of the six, the three plotted with points are the chaotic 
trajectory and the three plotted with lines are the non-chaotic trajectory. The "lower" two 
curves (starting with sibling distances near 10^^) arc integrated in extended precision with 
Taylor 1.4- The chaotic one fits an exponential curve with a Lyapunov time of about 12 
million years, while the non-chaotic one has the two trajectories separating approximately 
as t^-\ 
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Figure [8] plots two of these 21 systems. The value of lS.au jay differs between the two 
systems by one part in 10^. One of the systems appears chaotic, and the other does not, over 
a 200My timespan. The non-chaotic one has a semi-major axis of a;/ + 2 x 10~^, while the 
chaotic one has semi-major axis a;/ + 4 x 10~^. All other ICs in the two systems are identical. 
To ensure that the result is not integrator dependent, we have repeated the integrations with 
the Wisdom- Holman mapping included in Mercury 6.2 with timesteps of 8 and 4 days; and 
with NBI with a timestep of 4 days. As can be seen, all the integrations agree quite well 
with one another. Note that the Taylor I.4 integrations provide 3 extra digits of precision, 
and so the curves for Taylor I.4 are displaced about 3 orders of magnitude below the curves 
computed in double precision. Otherwise the shapes of the curves are virtually identical. We 
note that the chaotic one has a Lyapunov time of about 12 million years, while the regular 
one has the sibling trajectories separating from each other polynomially in time as t^'^. 



3.4. Accurate integrations over the age of the Solar System 



^The author has reported related results for integrations lasting 10^ years (IHayes. Wayne B. 

20071 ). The essential conclusion is the same, in that even after 10^ years, there remain some 



ICs (about 10%) that show no evidence of chaos, although some of the ICs appearing as 
regular over 200My develop exponential divergence later. 

Figure [9] displays the sibling divergence over 5 x 10^ years of the "canonical" IC used by 
DE405 (JED 2440400.5, June 28, 1969). As we can see, this IC shows little evidence of chaos 
for about the first 1.5Gy, and then develops slow exponential divergence with a Lyapunov 
time between about 200 and 400 million years. The individual planets each show similarly- 
shaped divergence curves (not shown), with the magnitude of divergence increasing with 
orbital radius. After 5 Gy, the uncertainty in Jupiter's position for this IC is less than 1 AU, 
while the uncertainty in Neptune's position is about 9 AU. Thus, there is a non-negligible 
chance that, if the Solar System lies close enough to the "canonical" IC of DE405, that we 
can know within about 10-15 degrees where each outer planet will be in its orbit when the 
Sun ends its main-sequence lifetime and becomes a red giant. Note that the levelling-off that 
starts at about the 4 Gy mark is not saturation (which occurs closer to 100 AU separation, 
while the separation at 5 Gy here is less than 10); the outer Solar System instead seems to 
be entering again into a period of polynomial (non-chaotic) divergence. 
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Sibling divergence plot of Outer Planets over 5 Gy 
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Fig. 9. — Sibling divergence over 5 Gy of the canonical IC of DE405, integrated using 
Taylor- 1.4. 



3.5. Percentage of initial conditions displaying chaos 

Observing the distance between sibling trajectories in Figures E] and at the 200My 
mark, we can reasonably simplify the distinction between chaotic and regular trajectories. 
By choosing a cutoff distance C and restricting our view to the various double-precision 
integrations, we can claim that siblings differing by less than C after 200My are regular (in the 
sense of having no observable Lyapunov exponent), while those differing by more than C are 
chaotic with a measurable positive Lyapunov exponent. Table |2] lists the number of systems 
that are chaotic by the above definition, as a function of timestep and cutoff. In addition 
to the 21-sample group of ICs displayed in Figures [5] and [6l we also drew a second sample 
group of 10 samples, spaced at 10- year intervals from 1900 to 1990. Since the testeph.f 
program included in DE405 provides only 7 digits of precision (corresponding to the accuracy 
to which the positions are known), taking ICs from DE405 at different times effectively takes 
ICs from different exact orbits, differing from each other by as much as one part in 10^. (As 
noted in a previous footnote, this method of sampling may not ideally represent an unbiased 
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sample from the observational error ball; instead, it is an unbiased sample from the set of 
7-digit-rounded initial conditions drawn from DE405.) We make several observations. First, 
the fraction of sampled systems that are chaotic by this simple definition is roughly about 
(70±10)%, and is relatively independent of both the timestep and the two sample groups (30- 
day vs. 10-year samples), although it of course increases as we decrease the cutoff. Recall that 
for a given initial condition, different stepsizes can give different results, so that which systems 
are chaotic changes as the timestep changes. However, here we measure only the number of 
chaotic systems as a function of timestep. If chaos were being "injected" into the system by 
the integrator, we would expect that the number of systems displaying chaos should increase 
with increasing timestep. However, this is not observed. A possible interpretation of this is 
that even a 400 day timestep reliably determines whether some system is chaotic, but the 
system is slightly differ ent for each timestep due to the IC perturbat ion introduced by the 



19921: 



Wisdona et al.lll996[) . Verifying this would 



symplectic integration ( Saha and Tremainel 

require us to implement "warmup" ( Saha and Tremainel 1993l or higher-order symplectic 
correctors than are included with Mercury 6.2 ( Chambera 1999 ). which we have not done. 
However, the fact that the fraction of systems displaying chaos is independent of timestep 
argues against the "chaos is injected by the integrator" hypothesis, at least for the timesteps 
used in this paper. 
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dt 


/lO 


C=l 
/21 


total 


% 


/lO 


C=0.1 
/21 


total 


% 


400 


7 


14 


21 


67.7 


9 


16 


25 


80.6 


200 


7 


11 


18 


58.1 


8 


15 


23 


74.2 


100 


6 


11 


17 


54.8 


8 


13 


21 


67.7 


050 


6 


14 


20 


64.5 


8 


15 


23 


74.2 


032 


7 


10 


17 


54.8 


9 


15 


24 


77.4 


016 


7 


14 


21 


67.7 


8 


16 


24 


77.4 


008 


8 


13 


21 


67.7 


8 


14 


22 


71.0 


004 


8 


13 


21 


67.7 


8 


18 


26 


83.9 


002 


8 


13 


21 


67.7 


8 


15 


23 


74.2 



Table 2: The percentage of initial conditions drawn from DE405 that admit chaos. Column 
labels: C is the cutoff described in the text; dt is the timestep; /lO means "out of the 10 
initial conditions drawn at 10-year intervals from 1900 to 1990"; /21 means "out of the 21 
initial conditions drawn at 30-day intervals starting at 1990" ; total is the sum of the previous 
two columns; % is the percentage of systems out of the 31 that display chaos according to 
the cutoff. 
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3.6. Chaos in the inner solar system is robust 



Varadi et al.l (120031 ) performed a 207My integration of the entire Solar System, including 
some non-Newtonian effects and a highly-tuned approximation to the effects of the Moon, 
and placed a lower bound of 30My on the Lyapunov time of the outer Solar System. However, 
they still saw chaos in the inner solar system. To test the robustness of chaos in the inner 
solar system, we performed several integrations of 8 planets (Mercury through Neptune) 
using the Wisdom-Holman mapping with timesteps of 8, 4, 2, 1, and 0.5 days. We treated 
the Earth-Moon system as a single body. To ensure that chaos in the outer solar system 
did not "infect" the inner solar system, we used only those DE405 ICs from Figures and 
E] for which the outer solar system was regular over 200My. We then integrated the system 
until chaos appeared. In all cases, even though the outer solar system was regular, the inner 
solar system displayed chaos over a short timescale such that information about the inner 
planetary positions was lost within about 20-50My. Thus, unlike the outer Solar System, 
we observe that chaos in the inner solar system is robust. 



4. Discussion 

We conclude that perturbations in the initial conditions (ICs) of the outer planets as 
small as one part in 10^ can change the behaviour from regular to chaotic, and back, when 
measured over a timespan of 200 million years. We believe this is the first demonstration 
of the "switch" from chaos to near-integrability with such a small perturbation of the ICs. 
Since our knowledge of the orbital positions of the outer planets is comparable to one part 
in 10^, it follows that, even if our simplistic physical model accounting only for Newtonian 
gravity were the correct model, it would be impossible at present to determine the Lyapunov 
time of the system of Jovian planets. Furthermore, it implies that an IC with 7 digits of 
precision (which is all an ephemeris can justifiably provide) can randomly lie on a chaotic 
or non-chaotic trajectory. Since our results converge in the limit of small timestep for the 
Wisdom-Holman mapping, and the converged results also agree with two very different 
high-accuracy integrations, and finally since the high-accuracy integrations in turn agree 
very well with quadruple-precision integrations, we believe that the results in this paper are 
substantially free of significant numerical artifacts. 



Guzzd ( I2OO5I ) corroborates the existence of a large web of 3-body resonances in the 



outer Solar System, and finds that their placement is consistent with Mu rray and Hol man's 



(1999) theory. Guzzo used his own fourth-order symplectic integrator (IGuzzd I2OOII ). and 
performed what appear to be reasonable convergence tests to verify the robustness of his 
main results. Thus, Murray and Holman's theory appears to explain the existence and 
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placement of 3-planet resonances. Furthermore, chaotic regions in Hamiltonian systems are 
usual ly densely packed with both chaotic and regular orbits (ILichtenberg and Lieberman 
19921 ). This paper corroborates the observation of densely packed regular and chaotic orbits, 
at a scale previously unexplored for the system of Jovian planets. 

As discussed in the text relating to Figure [8], we performed surveys across the semi- 
major axis au of Uranus in steps of 2 x lO"'^ AU for k = 6,7, 8. We found that, around the 
current best-estimate value of au, perturbations smaller than 2 x 10"^ AU had no effect on 
the existence of chaos. However, this does not imply that perturbations this small cannot 
have an effect; it simply means that the "border" between chaos and regularity is not within 
2 X 10~^ AU of the current best estimate of au. However, it is clear that the border between 
chaos and regularity is between the two systems depicted in Figure [HI which differ from each 
other in a^/ by 2 x 10^® AU. Although a survey across au for values between those two systems 
may not be very relevant from a physical standpoint, it might be very interesting from a 
dynamical systems and chaos perspective to probe the structure of the border between chaos 
and regularity. In p articular, it may be interesting to see if the border itself has some sort 
of fractal structure (jMandelbrotlll982l ). Suc h chaotic struc ture has already been observed in 
the circular restricted three-body problem (IMurisonl 119891 ). 



Newman et al. (2000) gave a compelling demonstration that the Wisdom + Holman 
symplectic mapping with too-large a timestep could introduce chaos into a near-integrable 
system by first showing that they could reproduce the chaos with a 400-day timestep, and 
then showing that the integration converged to being regular with a timestep of about 50 
days or less. Our results appear to hint that an even smaller timestep may be required: our 
curves depicting divergence-of-nearby-orbits did not fully converge until the timestep was 8 
days or less. On the other hand, the symplectic integrators produce solutions that effectively 
integrate a system w ith slightly perturbed ICs , although these pert urbations decrease with 
decreasing timestep (ISaha and Tremaind Il992l : IWisdom et al.l Il996l ) . Our Figure H] demon- 
strates that the behaviour does not always "switch" monotonically in timestep from chaotic 
to non-chaotic as observed by Newman et al. (2000); Table [2] also hints that the "amount" 
of chaos does not appear to increase with increasing timestep. An alternate interpretation 
is that even a 400-day timestep accurately integrates an orbit, but not the orbit that we 
chose. In particular, it may accurately integrate an orbit whose IC is perturbed slightly 
from the one we chose, and an appropriate correcti on may allow us to recover the correct or 



bit u sing an integration with much larger timestep (jSaha and Tremaindll992l : IWisdom et al. 



19961 ). However, the symplectic correctors in Mercury 6.2 are clearly not goo d enough to 
perform this recovery at a 400-day timstep; better correctors (jWisdonj 120061 ) or warmup 
(ISaha and Tremaind Il992l ) may be able to achieve this. 
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The convergence of our results at timesteps of 8 days or less, as well as the agree- 
ment with two different non-symplectic integrators (inclu ding the one used by New man 



Wisdom et al. 


Jl996[ 


are 


(Saha and Tremainel 


1992) 



our conclusions, although they might allow the same conclusions to be drawn using larger 
timesteps. 

With the exception of the results plotted in Figure [9], all of our simulations had a 
duration of 200 million years (My). All of them display an initial period of polynomial di- 
vergence, before the appearance of exponential divergence (if any). However, the duration of 
i nitial polynomial divergence differs greatly across systems, and has been observ ed by others 
( iLecar et al.ll200lh to last significantly longer that 200My; iGrazier et al.l (120051 ) observed it 
to last the entire duration of their 800My simulation. It would be interesting to create a table 
like Table [21 but including a "simulation duration" dimension as well. Certainly the evidence 
hints that more systems make the "switch" from polynomial to exponential divergence as 
the duration of the simulation increases. 

Our physical model is very simplistic, accounting only for Newtonian gravity between 
the Sun and Jovian planets. Although we ignore many physical effects which are known to 
effect the detailed motion of the planets (jLaskarlll999l : IVaradi et al.ll2003l ). it is unclear if such 
effects would substantially alter the chaotic nature of solutions. We at first believed that the 
largest such effect ignored was solar mass loss. Our first simulation s did not account for solar 
mass loss, which amounts to about one part in 10'' per million years (jLaskarlll999l : iNoerdlinger 
20051 ). Since we find that perturbations in position of that order can shift the system in- 
and-out of chaos, a naive analysis might lead one to suspect that solar mass loss might shift 
the planetary orbits in-and-out of resonance on a timesclase that is fast compared to the 
Lyapunov time, thus smoothing out the sibling divergence. We thus modified our model 
to include solar mass loss, but surprisingly it made absolutely no observable difference to 
any of the figures presented in this paper. To ensure that we did not make an error, we 
simulated systems with ever increasing mass loss until the Sun was losing 10% of its mass 
per lOOMy. We noted that the planetary orbital semi-major axes expanded significantly, as 
would be expected, but that the sibling divergences did not change until mass loss was at a 
rate of about 1% per lOOMy (1000 times greater than in reality). Thus, we conclude that 
solar mass loss also makes no difference to our results. 



-32- 



Acknowledgements 

The author thanks Scott Tremaine, Norm Murray, Matt Holman, PhiUp Sharp, and Bill 
Newman for helpful discussions and comments on the manuscript; Norm Murray and Philip 
Sharp for sending me (and explaining) their ICs; and Ferenc Varadi for giving me the source 
code to his most recent version of NBI. 



REFERENCES 

Barrio, R., F. Blesa, and M. Lara (2005). VSVO formulation of the Taylor method for 
the numerical solution of ODEs. Computers and Mathematics with Applications 50, 
93-111. 

Benettin, G. and A. Giorgilli (1994). On the Hamiltonian interpolation of near-to-the-identity 
symplectic mappings with application to symplectic integration algorithms. Journal 
of Stastistical Physics 7^(5/6), 1117-1143. 

Briggs, K. (1996). The doubledouble package: Software quadruple precision in C++. Un- 
published. 

Brouwer, D. (1937). On the accumulation of errors in numerical integration. Astron. J. 46, 
149-153. 

Chambers, J. E. (1999, April). A hybrid symplectic integrator that permits close encounters 
between massive bodies. Monthly Notices of the Royal Astronomical Society 304, 
793-799. 

Channell, P. J. and C. Scovel (1990). Symplectic integration of Hamiltonian systems. Non- 
linearity 3, 231-259. 

Grazier, K. R., W. I. Newman, J. M. Hyman, and P. W. Sharp (2005). Long simulations of 
the solar system: Brouwer's law and chaos. Preprint. 

Grazier, K. R., W. I. Newman, W. M. Kaula, and J. M. Hyman (1999, August). Dynamical 
Evolution of Planetesimals in the Outer Solar System. I. The Jupiter/Saturn Zone. 
Icarus 140, 341-352. 

Grazier, K. R., W. I. Newman, W. M. Kaula, F. Varadi, and J. M. Hyman (1995, May). 
An Exhaustive Search for Stable Orbits between the Outer Planets. Bulletin of the 
American Astronomical Society 27, 829. 



-33- 



Grazier, K. R., W. I. Newman, F. Varadi, D. J. Goldstein, and W. M. Kaula (1996, June). 
Integrators for Long-Term Solar System Dynamical Simulations. Bulletin of the Amer- 
ican Astronomical Society 28, 1181. 

Guzzo, M. (2001). Improved leap-frog symplectic integrators for orbits of small eccentricity 
in the perturbed Kepler problem. Celestial Mechanics and Dynamical Astronomy 80, 
63-80. 

Guzzo, M. (2005, March). The web of three-planet resonances in the outer Solar System. 
Icarus 174, 273-284. 

Hayes, Wayne B. (2007). Is the Outer Solar System Chaotic? Nature Physics. Accepted. 

Herbst, B. M. and M. J. Ablowitz (1989). Numerically induced chaos in the nonlinear 
schroedinger equation. Physical Review Letters 62, 2065-2068. 

Ito, T. and K. Tanikawa (2002, October). Long-term integrations and stability of planetary 
orbits in our Solar system. Monthly Notices of the Royal Astronomical Society 336, 
483-500. 

Jorba, A. and M. Zou (2005). A software package for the numerical integration of ODEs by 
means of high-order Taylor methods. Experimental Mathematics 14, 99-117. 

Kinoshita, H. and H. Nakai (1996). Long-Term Behavior of the Motion of Pluto over 5.5 
Billion Years. Earth Moon and Planets 72, 165-173. 

Laskar, J. (1989, March). A numerical experiment on the chaotic behaviour of the solar 
system. Nature 338, 237-238. 

Laskar, J. (1990, December). The chaotic motion of the solar system - A numerical estimate 
of the size of the chaotic zones. Icarus 88, 266-291. 

Laskar, J. (1994). Large-scale chaos in the Solar System. Astronomy and Astrophysics 287L, 
9-12. 

Laskar, J. (1996). Large Scale Chaos and Marginal Stability in the Solar System. Celestial 
Mechanics and Dynamical Astronomy 64, 115-162. 

Laskar, J. (1997). Large scale chaos and the spacing of the inner planets. Astronomy and 
Astrophysics 317, L75-L78. 

Laskar, J. (1999, August). The limits of Earth orbital calculations for geological time-scale 
use. Royal Society of London Philosophical Transactions Series A 357, 1735-1759. 



-34- 



Lecar, M., F. A. Franklin, M. J. Holman, and N. W. Murray (2001). Chaos in the Solar 
System. Annual Review of Astronomy and Astrophysics 39, 581-631. 

Lessnick, M. K. (1996). Ph. D. thesis, UCLA. 

Lichtenberg, A. J. and M. A. Lieberman (1992). Regular and Chaotic Dynamics. Springer. 

Lissauer, J. J. (1999). Chaotic motion in the solar system. Reviews of Modern Physics 71 (3), 
835-845. 

Mandelbrot, B. (1982). The Fractal Geometry of Nature. W. H. Freeman & Company. 

Morrison, L. V. and D. W. Evans (1998, November). Check on JPL DE405 using modern 
optical observations. A&AS 132, 381-386. 

Murison, M. A. (1989, December). The fractal dynamics of satellite capture in the circular 
restricted three-body problem. A J 98, 2346-2359. 

Murray, N. and M. Holman (1999, March). The Origin of Chaos in the Outer Solar System. 
Science 283, 1877-1881. 

Murray, N. and M. Holman (2001, April). The role of chaotic resonances in the Solar System. 
Nature 410, 773-779. 

Newman, W. I. and A. Y. Lee (2005, May). Symplectic Integration Methods and Chaos: 
Timestep Selection and Lyapunov Time. AAS/Division of Dynamical Astronomy 
Meeting 36. 

Newman, W. I., F. Varadi, A. Y. Lee, W. M. Kaula, K. R. Grazier, and J. M. Hyman (2000, 
May). Numerical Integration, Lyapunov Exponents and the Outer Solar System. 
Bulletin of the American Astronomical Society 32, 859. 

Noerdhnger, P. D. (2005). Solar mass loss, the astronomical unit, 

and the scale of the solar system. Preprint available at 

http : //home . comcast . net/ pdnoerd/SMassLoss . html. 

Rauch, K. P. and M. Holman (1999, February). Dynamical Chaos in the Wisdom-Holman 
Integrator: Origins and Solutions. AJ 117, 1087-1102. 

Reich, S. (1999). Backward error analysis for numerical integrators. SI AM J. Numer. 
Anal 36 {5), 1549-1570. 

Saha, P. and S. Tremaine (1992, October). Symplectic integrators for solar system dynamics. 
AJ 104, 1633-1640. 



-35- 



Sanz-Serna, J. M. (1992). Symplectic integrators for Hamiltonian problems: an overview. In 
A. Iserles (Ed.), Acta Numerica 1992, pp. 243-286. Cambridge University Press. 

Standish, E. M. (1998). JPL planetary and lunar ephemerides, 

DE405/LE405. JPL lOM 312.F-98-048 . August 26, 1998. Available at 
http : //ssd . jpl . nasa . gov/ iau-comm4/relateds . html. 

Sussman, G. J. and J. Wisdom (1988, July). Numerical evidence that the motion of Pluto 
is chaotic. Science 24 i, 433-437. 

Sussman, G. J. and J. Wisdom (1992). Chaotic evolution of the solar-system. Science 257, 
56-62. 

Varadi, P., B. Runnegar, and M. Ghil (2003, July). Successive Refinements in Long-Term 
Integrations of Planetary Orbits. The Astrophysical Journal 592, 620-630. 

Wisdom, J. (2006, April). Symplectic Correctors for Canonical Heliocentric n-Body Maps. 
AJ 131, 2294-2298. 

Wisdom, J. and M. Holman (1991, October). Symplectic maps for the n-body problem. The 
Astronomical Journal 102, 1528-1538. 

Wisdom, J., M. Holman, and J. Touma (1996). Symplectic Correctors. Fields Institute 
Communications, Vol. 10, p. 217 10, 217. 



This preprint was prepared with the A AS lAT^jX macros v5.2. 



